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Abstract 

We develop a numerical solver for the integral-differential equations, which de- 
scribe the radiative transfer of photon distribution in the frequency space with 
resonant scattering of Lya photons by hydrogen gas in the early universe. The 
time-dependent solutions of this equation is crucial to the estimation of the ef- 
fect of the Wouthuysen-Field (WF) coupling in relation to the 21 cm emission and 
absorption at the epoch of reionization. However, the time-dependent solutions of 
this equation have not yet been well performed. The resonant scattering leads to 
the photon distribution in the frequency space to be piecewise smooth containing 
sharp changes. The weighted essentially nonoscillatory (WENO) scheme is suitable 
to handle this problem, as this algorithm has been found to be highly stable and 
robust for solving Boltzmann equation. We test this numerical solver by 1.) the 
analytic solutions of the evolution of the photon distribution in rest background; 

2. ) the analytic solution in expanding background, but without resonant scattering; 

3. ) the formation of local Boltzmann distribution around the resonant frequency 
with the temperature to be the same as that of atom for recoil. We find that the 
evolution of the photon distribution due to resonant scattering with and without 
recoil generally undergoes three phases. First, the profile of the photon distribution 
is similar to the initial one. Second, an extremely flat plateau (without recoil) or 
local Boltzmann distribution (with recoil) form around the resonant frequency, and 
the width and height of the flat plateau or local Boltzmann distribution increase 
with time. Finally, the distribution around the resonant frequency is saturated when 
the photons from the source is balanced by the redshift of the expansion. This result 
indicates that the onset of the W-F coupling should not be determined by the third 
phase, but by the time scale of the second phase. We found that the time scale of 
the W-F coupling is equal to about a few hundreds of the mean free flight time 
of photons with resonant frequency, and it basically is independent of the Sobolev 
parameter if this parameter is much less than 1 . 
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1 Introduction 

It is generally believed that detecting redshifted 21 centimeter signals from 
early universe is one of the next frontiers in observational cosmology, because 
it would be able to provide the information of the first generation of light 
sources in the cosmic dark ages. Many studies have been done on the 21 cm 
emission and absorption from the halo of individual first stars (Chuzhoy et 
al. 2006, Cen 2006, Liu et al. 2007). A common conclusion of these works 
is that the configurations of the 21 cm emission and absorption regions is 
strongly time- dependent. The reason is simple. The necessary conditions of 
21 cm emission and absorption are that 1. the fraction of neutral hydrogen 
(HI) is still high; 2, Lya photons are available for the Wouthuysen-Field (W- 
F) coupling. The region of 21 cm emission around first stars is then a thin 
shell just outside the I-front (ionizat ion-front). Hence the 21 cm emission shell 
should move with a speed higher than the speed of the I-front Vf, which is 
rather high, even comparable to the speed of light. Therefore, the time scale 
of the formation and evolution of the regions of 21 cm signal can roughly be 
estimated by d/vf, d being the thickness of the 21 cm emission and absorption 
shell. This time scale is found to be of the order of 1 Myr and less (Cen 2006, 
Liu et al. 2007). 

In all of the above-mentioned works, the spin temperature is calculated with 
the assumption that the Wouthuysen-Field (W-F) coupling (Wouthuysen, 
1952; Field, 1958, 1959) is effective. That is, the color temperature of photons 
around Lya frequency is assumed to be the same as the kinetic temperature 
of hydrogen gas. The W-F coupling locks the color temperature of Lya pho- 
ton to the kinetic temperature of hydrogen gas, and then links the internal 
(spin-) degree of freedom with the kinetic temperature of gas. Obviously, the 
W-F coupling would be available if the time scale of the onset of the W-F 
coupling is less than the time scale of the formation and evolution of the 
21 cm emission and absorption shells. However, most calculations on the W- 
F coupling are based on the time-independent Fokker-Planck approximation 
(Chen & Miralda-Escude, 2004; Hirata, 2006; Furlanetto & Pritchard 2006; 
Chuzhoy & Shapiro 2006) of the radiative transfer with resonant scattering 
of Lya photons by hydrogen gas. The time-independent solution is acceptable 
only if the time-independent state is approached in a scale shorter than that 
of the evolution of the 21 cm signal regions. Unfortunately, this assumption is 
not obvious. Time-dependent solutions are necessary. 
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The W-F coupling is due to the resonant scattering of Lya photons by hydro- 
gen atoms. It is described by a Boltzmann-like integro-differential equation of 
radiative distribution in the phase space. The W-F coupling is sensitive to the 
evolution of the photon distribution in the frequency space. The analytical 
solutions of this equation with resr background has revealed that the fre- 
quency distribution of photons under the resonant scattering of Lya photons 
by gaseous hydrogen atoms without recoil has two features: 1. an extremely 
flat top of the photon profiles in the frequency space; 2. a sharp boundary 
of the flat top range(Field 1958). Currently, the numerical results are still far 
from precise to match these features (e.g. Meiksin, 2006). 

In this paper, we develop a numerical solver with the weighted essentially 
nonoscillatory (WENO) scheme. The WENO method has high order of ac- 
curacy and good convergence in capturing discontinuities as well as to be 
significantly superior over piecewise smooth solutions containing discontinu- 
ities (Shu 2003). WENO schemes have been widely used in applications. It 
is also effective in solving Boltzmann equations (Carrillo et al. 2003, 2006) 
and radiative transfer (Qiu et al. 2006, 2007, 2008). Therefore, one can expect 
that the integral-differential equations of resonant scattering can be properly 
handled numerically by the WENO scheme. 

The paper is organized as follows. Section 2 presents the basic equations of 
the resonant scattering of radiation. Section 3 gives the numerical solver of 
the WENO scheme. Section 4 presents the tests of the numerical solver. The 
time scale of the W-F coupling is briefly addressed in Section 5. A discussion 
and conclusion are given in Section 6. 



2 Basic equations 

2.1 Radiative transfer equations with resonant scattering 

Considering a spatially homogeneous and isotropically expanding infinite medium 
consisting of neutral hydrogen with temperature T, the kinetics of photons in 
the frequency space is described by the radiative transfer equation with reso- 
nant scattering (Hummer & Rybicki, 1992; Rybicki & Dell'antonio 1994) 

dJ ^ t) +2HJ( X ,t)- C -^ dJ ^ t) = 
Ot Vt ox 

-kc(f)(x)J(x, t) + kc [ TZ(x, x')J(x', t)dx' + C(t)<j>(x) (1) 
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where J is the specific intensity in terms of the photon number, H(t) = 
R(t)/R(t) is the Hubble parameter, R(t) is the cosmic factor, vt = {2UbT /m) 1 ' 2 
is the thermal velocity of hydrogen atom, the dimensionless frequency x is re- 
lated to the frequency v and the resonant frequency u by x = (u — z/ )/Az/ D , 
and Aud = vqVt/c is the Doppler broadening. The parameter k = x/Avd, 
and the intensity of the resonant absorption x is given by x = ' K e 2 nif 12/ m e c, 
where n\ the number density of neutral hydrogen HI at ground state, and 
/12 = 0.416 is the oscillation strength. The cross section at the line center is 

7TP 2 

do = fi2(Au D )-\ (2) 



In eq.flU), C(t) is the source of photons with the frequency distribution 4>(x), 
which is the Voigt function of the frequency profile with the center at Vq, i.e. 



-v 2 



{X) = ^J dy (x-y) 2 + a 2 (3) 



-oc 



where a is the ratio of the natural to the Doppler broadening. We have a = 
A2i/{8irAb>D), and A21 = 6.25 x 10 8 Hz is the Einstein spontaneous emission 

coefficient. 0(x) is normalized with J ^x')dx' = 1. When a _ 0, we have 

pure Doppler broadening as 

<> D (x) = -^e~ x2 . (4) 

'7T 



The redistribution function TZ(x, x') gives the probability of a photon absorbed 
at the frequency x', and re-emitted at the frequency x. It depends on the details 
of the scattering (Henyey 1941; Hummer 1962; Hummer, 1969). If we consider 
coherent scattering without recoil, it is 

n(x,x')= (5) 



1 



00 



^3/2 

\x-x'\/2 



e- 2 



tan — tan 1 



du 



where = m i„(^0 and _ = ^,,0. 0bv i0U8ly , jll^W ^ 

— OO 

4>(x). In the case of a = 0, i.e. considering only the Doppler broadening, eq.Q 
becomes 

lZ(x,x') = -erfc[max(|x|, \x'\)]. (6) 
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Considering the recoil of atoms, the redistribution function for the Doppler 
broadening is (Field, 1959, Basko, 1981) 

TZ{x,x') = ^e 26x ' +b2 erfc[max(|x + b\, \x' + 6|)], (7) 



where parameter b = hu /mvTC = 2.5 x 10 4 (10 4 /T) 1//2 . 



2.2 Re-scaling the equations 



We use the new time variable r defined as r = cnia^t, which is in units 
of the mean free flight time of photons at resonant frequency. The number 
density of neutral hydrogen atoms n x = fm n H, with /hi being the fraction of 
neutral hydrogen. For the concordance ACDM model, we have hh = 1.88 x 
10~ 7 (Qbh 2 /0.02)(l + z) 3 cm~ 3 . The factor 0.75 is from hydrogen abundance. 
Therefore, 

, T v l/a , 1Q v a / .022\ 
^0.0Mrf n l( w ) (_) (^) yrs (8) 

We re-scale the eq.((T]) by the following new variables 

J'(x, t) = [R(t)/R } 2 J, C\t) = [R(t)/R ] 2 C. (9) 



Thus, eq.([T]) becomes 



dJ'(x,r) , 

= -<f>(x)J (X,T) 

r 8 7' 

+ J K(x, x')J'(x', r)dx' + 7 — + C\t)4>{x). (10) 

We will use J for J 1 and C for C in the equations below. The parameter 
7 is the so-called Sobolev parameter, 7 = (H/vxh) = (8nH/3Ai2\ 3 ni) = 
(iYm e z/o/7re 2 ni/i2), where A is the wavelength for Lya transition. 7 is simply 
related to the Gunn- Peterson optical depth tqp by 

7 - = tgp = 4.9 x f^V' 2 (^f ■ (11) 



Qm/ V 0.022 } V 10 



The redshift evolution of /hi is dependent on the reionization models. Before 
reionization /hi — 1; after reionization /hi — 10~ 5 in average. Therefore, the 
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parameter 7 has to be in the range from 1 (z < 7) to KT 7 (z > 10). For static 
background, 7 = 0. 



3 Numerical solver: the WENO scheme 

3.1 Computational domain and computational mesh 

The computational domain in the case of static background is x G [—6,6]. 
The initial condition is J(x, 0) = 0. The boundary condition is 



In the case of expanding background, i.e. 7 7^ 0, the computational domain is 
bigger than x G [—6, 6] depending on the value of the Sobolev parameter 7. 
The domain {xi e f t ,x r i g ht) is chosen such that for the particular value of 7 we 
have 



For example, the domain is taken to be x G [—100, 6] for the case of 7 = 10 -3 . 
Also for different values of 7 the solutions reach saturation at different time. 
For example, we see in our numerical results that the solution J(V, t) of eq. 
(Q reaches saturation at time r = 100 for 7 = 10 _1 , and reaches saturation 
at time r = 10 4 for 7 = 10~ 3 . The computational domain is discretized into a 
uniform mesh in the x direction, 



where Ax = (x rig ht—xieft)/N x , is the mesh size. We also denote J™ = J(x«, r n ), 
the approximate solution values at (xj,r™). 



3.2 Algorithm of the spatial derivative 

To calculate we use the fifth order WENO method (Jiang k, Shu, 1996). 
That is, 



J(x, t) = 0, at \x\ — 6. 



(12) 



J(xi e ft,T) W 0, J{x r ight,l~) ~ 0. 



(13) 



Xi = xi e ft + iAx, 



i = 0,l,2,.-.,JV s 



dJ{Xi,T n ) 



1 

Ax 



{hj+1/2 ~ hj-1/2) 



(14) 



dx 
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where the numerical flux is obtained by the procedure given below. We 

use the upwind flux in the fifth order WENO approximation because the wind 
direction is fixed (negative). First we denote 

hi = J( Xi ,T n ), i = -2, -1,..., N x + 3 (15) 



where n is fixed. The numerical flux from the WENO procedure is obtained 
by 

k+l/2 = + ^2^+1/2 + ^3^1/2, ( 16 ) 



where h\™y 2 are the three third order fluxes on three different stencils given 
by 



1 5 X 

f (3) _ 11 , 7 1 

"j+1/2 — ~Q i+1 ~ Q i+2 3 *+ 3 ' 

and the nonlinear weights u m are given by, 

^--3—, ^"(TTaP ( } 



where e is a parameter to avoid the denominator to become zero and is taken 
as e = 1CT 8 . The linear weights 7? are given by 

71 =10' 72 = 5' 73 =10' (18) 



and the smoothness indicators fli are given by 



13 1 
Pi = 3^(^-1 - 2h i + ^+i) 2 + 4(^-1 - 4/i * + 3^+i) 2 , 

13 1 

#2 = — (hi - 2h i+1 + h i+2 ) 2 + -(hi - h i+2 ) 2 , 

13 1 

P3 = Y2^ hi+1 ~ 2hi+2 + hi+ ^ + 4^ 3/li+1 ~ 4/li+2 + ^+3) 2 - 
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3.3 High order numerical integration 



The integration of the resonance scattering term is calculated by a fifth order 
quadrature formula (Shen et al. 2007) 



f f(x)dx = Axjrwjfixj) + 0(Ax 5 ), (19) 

3=0 



Xleft 



where the weights are defined as, 



251 299 211 739 

w = ^7T' w l = 7T77T> w 2 = 7T77T> w 3 



720' 240' 240' 720' 

739 211 299 251 

W ^- 3 = 720' W ^- 2 = 240' WN '~ 1 = 240> WN * = 720' 



and Wj = 1 otherwise. Notice that this numerical integration is very costly 
and uses most part of the CPU time. For the equation with recoil and the 
redistribution function with a = 0, we have used a grouping of the numerical 
integration operations at different x locations so that the computational cost 
can be reduced to order O(N) rather than 0(N 2 ), where N is the number of 
grid points in x, without changing mathematically the algorithm and its accu- 
racy. The O(N) algorithm for numerical integration is given in the appendix 
which further highlights the speed and accuracy of the numerical algorithm 
proposed in this paper. Unfortunately, this grouping technique does not work 
for the case a ^ 0, hence the CPU cost for the case with a ^ is much larger. 



3.4 Time evolution 



To evolve in time, we use the third-order TVD Runge-Kutta time discretiza- 
tion (Shu & Osher, 1988). For systems of ODEs u t = L(u), the third order 
Runge-Kutta method is 

H (1) =« n + ArL(«",T n ), 

«(2) = ^ u n + ^(« (1) + AtL(u {1 \ r n + Ar)), 

u (3) = \ u n + 2 (u (2) + ArL ( u (2) > T n + 1 Ar)) _ 
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4 Tests of the WENO solver 



4-1 Static background 



We first test the WENO solver with two analytical solutions of eq. (TIP]) with 
static background, H = 7 = and Doppler broadening a = (Field, 1958). 
The first analytical solution is for the initial radiative field J(x, r) = and 
the constant source C(r) = 1. It is 



J{x, t) = tt- 1/2 [1 - exp(-re-* 2 )] (20) 
+ / e w2 [l - (1 + Te~ w2 ) exp{-Te- w2 )}eii{w)alw. 



The second analytical solution of eq. ffTUj) is also for H = 7 = 0, a = 0, but the 
source C = 0, while the initial radiative field is 

J(y,0) = 7r- 1 / 2 e^ 2 . (21) 



The solution is 



00 

J(x, t) = Ti-^e-* 2 exp{-Te~ x2 ) +t J e"^ exp(-r e - w2 )erf (w)dw. (22) 



The analytical solutions (1201) and (j22p are shown in Figures [T] and [2] respec- 
tively. We also plot the numerical solutions given by our algorithm in same 
figures. The numerical results show very small deviation from the analytical 
solutions. 

A common feature of Figures 1 and 2 is that the originally Doppler peak at the 
center of the frequency profile gradually becomes a flat plateau. The width 
of the plateau increases with time. This is because the resonant scattering 
makes a non-uniform distribution in the frequency space to a uniform one. It 
is similar to that in the physical space, diffusion generally leads to an evolution 
from non-uniform distribution to a uniform one. The height of the plateau of 
Figure 1 is increasing with time, while in Figure 2 it is decreasing, because for 
the case C = 1, the number of photons increases, while for the case of C = 0, 
the total number of photons is conserved. 
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Fig. 1. Static solutions (7 = 0) for pure Doppler redistribution (a = 0) of eq. lfTUj) . in 
which C = 1 and J(x,0) = 0. The analytical solutions are shown by dashed lines, 
while the numerical results are shown by solid lines. 
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Fig. 2. Static solutions (7 = 0) for pure Doppler redistribution (a = 0) of eq, (fTUj) . 
in which C = and J(x, 0) = 7r~ 1//2 exp(— x 2 ). The analytical solutions are shown 
by dashed line, while the numerical results are shown by solid lines. 



4-2 Expanding background 



The second test is given by considering expanding background. Without ab- 
sorption and scattering, eq.(10) becomes 

« =7 | +0(tWx) . (23) 
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Fig. 3. Solutions of J(x, r) of eq. ([23j) with initial condition J(x, 0) = 0. The param- 
eter 7 is taken to be 10 _1 (left panel) and 10 -3 (right). The analytical solutions 
eq.(24) are shown by dashed lines 
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Fig. 4. Solutions of J(x, r) of eq. (|10p with initial condition J(x, 0) = 0. The param- 
eter 7 is taken to be 1CP 1 (left) and 10 -3 (right). 



If C is r- independent, the analytic solution of eq.(23) is (Rybicki & Dell'antonio, 
1994) 



J(x, t) = J(x + jt, 0) + Cj^l^x) - $(x + 7r)] 



(24) 



where 



$(x) = / <p(x')dx' 



(25) 
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The meaning of the solution eq.(24) is simple. It just shows that the redshift 
of photons in frequency space is described by 



— x = 77-. 



(26) 



We take <p(x) = (l/^)e~ x , and assuming the initial field J(x, 0) = 0. Figure 
3 shows that the numerical results obtained by using the WENO algorithm 
agree well with the analytic solution of eq.(24). 

Figure 4 is the solutions of eq.(10) with parameter 7 = 10 _1 and 10 -3 . In 
these cases the analytical solutions are not available. We can see from Figure 
3 and 4 that the intensity J stops to increase and approaches a saturated value 
when r is large. This happens when the number of photons produced from 
the sources is equal to that of the redshifted photons. 

4-3 The effect of recoil of atom 



Fig. 5. Solutions of J(x, r) of eq. (Tlu]) with redistribution function eq.(7), and 7 = 0. 

The third test is by using the redistribution function with recoil Eq.(7). Fig- 
ure 5 gives the solutions of eq.(10) with parameter b = 0.03, and the initial 
condition and source term are taken to be the same as Figure 1. The results 
of Figure 5 actually is similar to Figure 1, but only the center flat plateau of 
Figure 1 now is replaced by a sloping plateau. The latter is a local Boltzmann 
distribution around the renasont scattering. 

Field (1959) has shown that, once the solution J(x,r) of eq. ljTOl) with 6 = 
has a flat plateau in the central region, the effect of recoil is to make the flat 
plateau to be a Boltzmann-like distribution, i.e. if the solution with no recoil 
redistribution eq.([6]) shows J(x,r) ~ J(0,r), within \x\ < Xf, the solution of 



T= 10000 




X 
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eq.(10) with redistribution eg. (123!) will be 

J(0, r)e' 2bx = J(0, r)e- h{u -" o)/kT , \x\ < x f . (27) 

One can test this property with F{r) defined as 

F(r) = logJ(0,r)-log(l,r) (28) 

If the local Boltzmann distribution within \x\ < 1 is realized at time r, F(t) 
should be equal to 2b [eq.(27)]. In Figure 6, we present the relation of F(r) 
vs. r. At t = 0, the Gaussian sources [eq.(21)] yields F(0) = 1. Then, F(r) 
approaches to 2b at r > 10 2 . These results show that the WENO solver is 
reliable and effective to study the time- dependence of J described by the 
resonant scattering equation. 

Figure 6 also shows the relation F(t) vs. r for solutions b = 0.03 (left panel) 
and b = (right panel). The two curves actually are similar. When the right 
curve approaches to F(r) — > at time r, the left curve, at the same time, 
approaches to F(r) — > 2b. Therefore, it would be reasonable to time-scale of 
the formation of local Boltzmann distribution by the formation of the flat 
plateau. 




Fig. 6. F(t) as function of r for solutions shown in Figure 5 (left) and Figure 1 
(right) 

5 The Wouthuysen-Field Coupling 

We now estimate the time scale, twf, of the onset of the W-F coupling. As 
mentioned in §4, in the first phase, J(x, r) keeps the initial profile (f>(x) of 
the photons from source C. In the second phase, the profile is no longer the 
initial one. A flat plateau (without recoil) or local Boltzmann distribution 
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(with recoil) form the central part (\x\ ~ 0) due to the resonant scattering. 
The width of the flat plateau increases with the time r. Finally, in the third 
phase, the injection of Lya photons from the source C is balanced by the 
redshift, the height of the flat plateau or the local Boltzmann of J(x, r) will 
stop increasing, and reach to a saturated value. The width of the flat or sloping 
plateau on the red side continuously increases. 

As mentioned in Section 1, most works on the effect of the W-F coupling of 
21 cm signal are based on the static solution of the Fokker-Planck approx- 
imation of eq.(10) (Chen & Miralda-Escude, 2004; Hirata, 2006; Furlanetto 
& Pritchard 2006; Chuzhoy & Shapiro 2006). Obviously, this solution corre- 
sponds to the third phase. Thus, the time scale of the W-F coupling is of the 
order of rm- Therefore, this approximation would be reasonable if the time 
scale of the formation and evolution of the 21 cm signal regions is larger than 
Tin- 

From Figures H] and [BJ one can see that the time scale of the onset of the 
third phase, Tin, has to be of the order of 10 2 , and 10 4 for 7 = 10 _1 , and 10~ 3 
respectively, i. e. r in is roughly equal to a few to ten of Gunn-Peterson optical 
depth. On the other hand, we can see from eq. ffTTl) that 7 = 10 _1 , and 10~ 3 
correspond to /hi = 10~ 4 , and 10~ 2 for the source at redshift 1 + z = 10. Thus, 
for all cases, eq.® yields that the time scale tm to be of the order of 1 Myr. 
The static solution of the Fokker-Planck equation would not be valuable for 
the 21 cm problem, if the time scale of the evolution of the 21 cm region is 
equal to or less than 1 Myr. 

Actually the time-independent solution of eq. ffTOl would not be necessary 
for the 21 cm signal. The relative occupation of the two hyperfme-structure 
components of hydrogen ground state depends only upon the shape of the 
spectrum near the Ly-alpha frequency, regardless whether the solution is time- 
independent, or saturated. What we need for the W-F coupling is only the 
frequency distribution to show a local Boltzmann-like distribution J(x, r) oc 
exp[— h(u — Uo)/kT] around vq, where T is the kinetic temperature of hydrogen 
gas. The formation of the Boltzmann-like distribution is irrelevant to the time 
scale Tm, but depends on the onset of the second phase of the J(x, r) evolution. 

Therefore, twf should be estimated by the time scale of the onset of the second 
phase, m, n °t by Tm. According to the property shown by Figure 6, tu can be 
found by the formation of a flat plateau around the resonant frequency. For 
small 7, Tu is much less than rjn- To show this point, we calculate the solution 
of eq.(10) with 7 = 10~ 5 and 10~ 7 on static background. The results are given 
in Figure 7. It shows that that a small flat plateau has already formed at 
r = 100. On the other hand, Tm will be as large as r ~ 10 6 for 7 = 10 -5 . 

In Figure El we show the function F(t) for the parameters (a = 0, 7 = 10~ 5 ), 
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Fig. 7. Solutions of J(x,t) of eq.(10) with 7 = 10 5 (left panel) and 10 7 (right 
panel). 





Fig. 8. F{t) as function of r for solutions of a 
10" 7 (right). 



and 7 = 10~ 5 (left) and 



and (a = 0, 7 = 10~ 7 ). It is interesting to observe that the two sets of F(t) 
are almost identical. That is, the time scale twf is actually independent of 
the parameter 7. For both cases of the parameter 7, twf is equal to about a 
few of 10 2 . In other words, a flat plateau with proper width can form with 
a few tens or hundreds of the mean free flight time of the resonant photons, 
regardless the expansion of the universe. 



The reason of the 7-independence of Tn or t^f can be directly seen from 
eq.(10). This equation describes two kinetic processes of approaching statisti- 
cally steady state. The first is the resonant scattering, which leads to the for- 
mation of a flat plateau or local Boltzmann distribution around the resonant 
frequency. The second is the redshift of photons, which leads to a statistical 
equilibrium between the injected photons and redshifted photons. The steady 
state due to resonant scattering corresponds to the second phase, while the 
steady state due to redshift corresponds to the third phase. Similar to various 
statistical equilibrium or steady state maintained by collision or scattering, 
the steady state of resonant scattering can be realized via a few tens or hun- 



15 



dreds of the resonant scattering. On the other hand, 7 is a ratio between 
the time scales of the resonant scattering and the expansion of the universe. 
Therefore, when 7 is much less than 1, the time scale of the formation of 
the flat plateau and local Boltzmann distribution are much shorter than the 
expansion of the universe. This point can also be seen with Figure [6l which 
shows that for 7 = 10~ 5 and 10~ 7 the evolutions of J(x, r) actually are the 
same till r = 10 4 . In short, the W-F coupling will take place once a statistical 
equilibrium localized around the resonant frequency is realized. 



6 Concluding remarks 

To study the time scale of the W-F coupling, we need a better solution of the 
integro-differential equation of the resonant scattering of Lya photons. Espe- 
cially, we need to know the evolution of photon frequency distribution around 
the resonant scattering. That is, the algorithm should be able to handle the 
extremely flat distribution and its sharp boundary. These features can prop- 
erly be captured by the WENO scheme, which has high order of accuracy and 
good convergence in capturing discontinuities as well as to be significantly su- 
perior over piecewise smooth solutions containing discontinuities. The WENO 
algorithm can be used for resonant scattering with and without the recoil of 
atoms. This algorithm is reliable as it passed all the tests. 

The evolution of the photon distribution in the frequency space generally un- 
dergoes three phases. In the first phase, the profile of the photon distribution 
is similar to the initial one. In the second phase, an extremely flat plateau 
(without recoil) or Boltzmann distribution (with recoil) formed around the 
resonant frequency. The width and height of the flat plateau or local Boltz- 
mann distribution increase with time. Finally, in the third phase, the photons 
injected from the source is balanced by the redshift of the expansion, and the 
evolution of the photon distribution is stable. The first phase is very short. 
The second phase will be onset after a few tens or hundreds of photons scat- 
tering by atoms. On the other hand, the onset of the third phase is mainly 
dependent on the Gunn-Peterson optical depth, which is large at early uni- 
verse. Consequently, the onset of the third phase is much later than the second 
one. 

Usually, the W-F coupling is described by time-independent solutions of the 
Fokker-Planck approximation of the integro-differential equation of the reso- 
nant scattering. With the WENO solutions, we show that the time-independent 
solutions would not be available for the 21 cm signal of the first generation 
stars, if the life time of the evolution of the 21 cm region is equal to or less 
than 1 Myr. However, the time scale of the onset of the W-F coupling actually 
is irrelevant with the time-independent solutions. The W-F coupling will take 
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place once a statistical equilibrium is locally realized in the frequency space 
around the resonant frequency. This time scale is always a few tens or hun- 
dreds of the mean free flight time of the resonant photons, and is generally 
independent of the expansion of the universe when the Gunn-Peterson opti- 
cal depth is large, or the Sobolev parameter is much less than 1. More detail 
results of the time-dependence o the W-F coupling is reported in Roy et al 
(2009). 
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A Numerical integration: an O(N) algorithm 

I(x) = ^J e 2te ' +b2 erfc(max(|x + 6|, \x' + b\))J{x',t)dx', (A.l) 

with lZ(x,x') as in eq.©. To evaluate I m = I(x m ), Vm = —N t , ■ ■ ■ ,N r , we 
apply the rectangular rule, which is spectrally accurate for smooth functions 
vanishing at boundaries, 

<** right 

I m = \ J erfc(max(|z m + 6|,|a/ + b\))e 2bx ' +b2 J(x',t)dx' (A.2) 

1 N r 

^-Ax £ erfc(max(|x m + 6|,|x, + 6|))e 2 ^ +62 J(x i ,t). (A.3) 

2 i=-N, 

Notice that this summation algorithm is very costly as it takes 0(N) opera- 
tions per m, therefore the total procedure has 0(N 2 ) operations overall. We 
use a grouping technique, described below, so that the overall computational 
cost can be reduced to O(N), without changing mathematically the algorithm 
and its accuracy. 

The proposed scheme with order N computational effort is the following. Let 
Nb = floor (■£-) and N 2 b = floor(j^-). The integration algorithm is designed 
for two cases: m > —Nb and m < —Nb. 

In the case of m > —Nb or equivalently x m + b > 0: 
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i -m-N 2b -l 

I m = -Ax( E erfc d^ + b\)e 2b ^ +b2 J(x t , t) 

1 i=-N, 

m 

+erfc(|a; m + 6|) E e 2hx ^ J(x h t) 

i=—m-N 2b 

N r 

+ J2 eric(\x t + b\)e 2bx > +b2 J(x t ,t)) (A.4) 

i=m+l 

= ^Ax(/ lim + erfc(|x m + b\)I 2 , m + h,m) (A. 5) 

(1) Evaluate h-N b , h-N b and h-N b respectively as 

N b -N 2b -1 

h,-N b = E erfc(|^ + &|)e 2te<+62 J(^,i), (A.6) 

i=-N t 
-N b 

h,-N b = E e 2b ^ +b2 J(x u t), (A.7) 

i=N b -N 2b 
N r 

h,-N b = E erfcd^ + feDe 2 ^^ 2 ^^), (A.8) 

i=-N b +l 

which leads to O(N) cost. 

(2) DO m = —N b + 1, N r Evaluate J 2im , respectively by 

h, m = h,m-i -erfc(|x_ m _^ 26 + b\)e 2bx — ^ J{x^ N2b ,t) (A.9) 
h, m = h,m-i + e 2hXm ^ J(x m , t) + e 2bx — + # J(x_ m _ N2b , t) (A.10) 
/ 3) m = / 3 , m -i -erfc(|a; m + 6|)e 2 ^ +62 J(a: m ,t) (A.ll) 
ENDDO 

To be consistent with the indeces, if Ni — N 2 b < N r then, we will set 
h,m = 0, for m = Ni — N 2 b, N r . The algorithm leads to 0(1) cost per m, 
therefore 0{N) computation overall. 

In the case of m < —N b , or equivalently x m + b < 0: 



i m— 1 

/ m = -Ax( E erfc(|^ + 6|)e 2 ^ +b2 J(^,t) 

-m— JV 2i) -l 

+erfc(|x m + 6|) E e 2bx * +b2 J( Xl ,t) 

i=m 

N r 

+ E erfc(|^ + 6|)e 2 ^ +62 J(^,0) (A.12) 

i=-m-No h 
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= ^Ax(I 1>m + erfc(|x m + b\)I 2%m + h, m ) (A. 13) 

(1) Evaluate h _jv 6 -i, ^2,-^-1 and i" 3 ,-JV 6 -i as 

-N b -2 

Ii,-n„-i= E erfc(|x l + &|)e 2 ^ +fe2 J(x l ,t), (A.14) 
V"*-i = E e^V^t), (A.15) 

i=-JV 6 -l 

V^-i= E erfc(|^ + 6|)e 2b3:i+62 J(^,0, (A.16) 

i=JVi,+l-JV 2 6 

which leads to 0(N) cost. 

(2) DO m = -N b - 2, - jV, 

Evaluate i"i, m , J 2jm , J 3jTn respectively by 

= /i,m+i - erfc(|x m + 6|)e 26 ^ +62 J{x m , t) (A.17) 
/ 2 , m = h, m+ i + e 2bx ™ +b2 J(x m ,t) + e 2bx —^ +b2 J(x_ m _ N2b ^,t)(A.18) 
h, m = h, m+ i-eTi C {\x_ m „ N2h ^ + b\)e 2bx -™- N 2 b -^ 
ENDDO 

To be consistent with the indeces, if Ni — N 2 b > N r , we will set I^ m = 0, 
for m = —N 2 b — N r — 1, —N t . Again, the algorithm leads to 0(1) cost per 
m, therefore O(N) computation overall. 
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